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suggested  for  modeling  time  series  data  that  contain  episodes,  j,.e.  large, 
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rate  parameter.  The  method  is  illustrated  using  hormone  levels  data 


1.  INTRODUCTION 


Let  (Xj_ ,  t*0,  +1,  +2,...)  be  a  discrete  time  series,  l.e.  the  xt's 
are  a  collection  of  random  variables  with  finite  second  moments.  In 
analyzing  time  series  data  one  usually  assumes  that  X  Is  covariance  stationary, 
l.e.  the  mean  value  function  mt  *»  E[Xt]  is  a  constant  and  that  the  covariance 
function  K(s,t)  «  Cov[X  ,X  ]  is  only  a  function  of  |s-t|.  Unfortunately, 

S  L 

such  an  assumption  rules  out  a  wide  class  of  interesting  series.  For 
example.  Figure  1  is  a  graph  of  levels  of  Luteinizing  hormone  (LH)  in  a 
cow  measured  at  ten  minute  intervals  for  24  hours.  It  is  difficult  to 
imagine  that  this  data  comes  from  a  covariance  stationary  time  series.  Even 
if  it  were  the  usual  methods  of  time  series  analysis  (correlation  and  spectral 
analysis)  would  not  address  the  questions  the  scientist  would  be  interested 
in,  namely  the  intensity  of  episodic  occurrence,  the  average  heighth  of 
episodes,  and  the  rate  at  which  the  hormone  level  returns  to  its  baseline 
value. 

The  features  that  distinguish  this  data  (and  other  data  in  a  wide 
variety  of  areas  such  as  river  flows,  pollution  levels,  or  chemical  manu¬ 
facturing  processes)  are  the  random  occurrence  of  large  surges  followed  by 
an  exponential  return  to  the  low  relatively  constant  baseline  values.  A 
model  which  reflects  such  behavior  was  suggested  by  Parzen  (1962b).  We 
assume  that  the  process  started  at  some  time  zero  and  let  the  number  of 
episodes,  N(t),  from  time  0  to  time  t,  t  i  0,  be  a  Poisson  process  with 
intensity  parameter  X,  i ,e .  for  any  n  1  and  any  times  0  t^  <  t2  K  •••  <tn* 
N(tn)  -  N(tn_^),...,  N(t2>  -  N(t^)  are  independent,  Poisson  distributed 
random  variables  with  means  X(t  -t  ,),...,  X(t,  -  t.,).  Then  the  value 
Xt  of  the  process  at  time  t  is  assumed  to  be  of  the  form 
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d.i) 


where  the  y  '  s  are  lid  random  variables  with  mean  p  and  variance  a2,  while 
*'m  y  y 

is  the  time  of  occurrence  of  the  mth  episode.  Thus  the  time  series  consists 
of  the  sum  of  a  random  number  of  exponentially  decaying  spikes  where  the 
number  of  such  spikes  is  governed  by  X,  their  heighth  by  p^,  and  the  rate 
of  decay  by  y.  The  aim  of  this  paper  is  to  provide  good  estimators  of  these 
parameters.  In  section  two  we  derive  the  properties  of  the  model  (1.1)  and 
in  section  three  we  find  consistent,  asymptotically  normal  estimators  of  X, 
liy,  and  y.  In  section  three  we  will  assume  that  y  is  large  enough  relative 
to  X  that  one  spike  decays  quickly  enough  before  the  occurrence  of  the  next  so 
that  they  do  not  overlap.  This  restriction  seems  reasonable  for  the  type  of 
data  sets  we  are  interested  in.  We  also  will  assume  that  the  baseline  value 
is  essentially  constant.  In  many  episodic  time  series  the  variation  in  the 
baseline  is  of  importance  also.  Any  method  of  analyzing  such  series  would 
have  to  be  iterative  in  nature,  first  considering  the  baseline  and  then  the 
episodes.  In  this  paper  we  consider  only  the  episodes  and  leave  the  general 
problem  to  future  research. 

2.  PROPERTIES  OF  THE  MODEL 
From  (1.1)  we  can  write 


Xt+1  S  ym  6  " 

m»»l 


-  T  y  T1’ 


tdfI 


m 


m*N(t)+l 


“  pXt  +  Gt+1 
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_l/v 

where  p  -  e  and  the  e^'s  are  iid  random  variables  since  they  are  a 
function  of  N(*)  for  the  nonoverlapping  time  intervals  (0,1],  (1,2],  etc. 
Further,  et+^  *s  independent  of  Xfl  for  s  <  t  +  1.  Thus  X  is  in  the  form 
of  an  autoregressive  process  of  order  one  but  as  we  shall  see  It  Is  not 
covariance  stationary.  We  note  that  for  0  <  y  <  ",  we  have  0  <  p  <1. 

Now  for  a  poisson  process  with  intensity  X  we  have  that  the  arrival 
times  in  the  interval  [s,t]  given  that  N(t)-N(s)  =  K  are  iid 

1  K 


uniformly  distributed  on  the  interyal  [s,t].  Thus 


E[£t1  EN(t)-N(t-l)  ]  E  . 


m=N (t-l)+l 


-(t-x  )/y 

ym  m  |n (t)-N(t-l) 
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*  EN(t)-N(t-l)[K  E<^>  E»Ct-l,t)<e"<t‘T“>/T>] 
■  EN(t)-N(t-l)  [K\ 


=  Xuy  y(l-e  1^Y)  , 


where  Eu(t_i  denotes  expectation  with  respect  to  the  Uniform  (t-l,t)  pdf. 


Further, 


E^|N(t)-H(C-l)  -  kJ  -  K  E(y2)Eu(t.l  t) 

K(K-l)  t«ym)l2[Eu(t.lit)(e-<t'T)/1']2 


K(u2+<>2)  y  (1-e  2/y)  +  K(K-l)  u2  y2(l-e~1^Y)2  , 
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and  so 


E[e2] 


X  (a2  +  y2)  1.  (i_e~2/Y)  +  y2  X2 
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Var[et]  -  \6±  (1-e  T) 


where  6  ■  y2  +  o2  . 

y  y 


Now  to  find  the  moments  of  Xfc  we  write 
Xt  -  pXt-l  +  £t 


P(pXt_2  +  Et_i)  +  et 


t  "v1  j 

>\ +  l  p  Vi  • 

i-0  c  2 


Thus 


t-1 


E[X  ]  =  pt  E[X  ]  +  Xy  Y(l-e~1/Y)  £  pj 

c  °  y  j-o 

■=  p*"EtX  ]  +  Xy  y(l-e"t/Y) 
o  y 


Also, 


t-1 


t+s-1 


“"'VW  ‘  Cov[pt  X„  +  pJet-j’pt+8Xo  +  ^Vs-J1 


t-1 


p2t+sVar[X  ]  +  Var[e  ]  pS  £  p2^ 

°  c  j-o 


-  p2t+S  Var [X  ]  +  X6*  (l-e"2/Y)  PSU-P2t)/(l-P2) 
o  z 


p2t+s  Var[X  ]  +  XS*  p8(l-p2t)  , 


'which  gives 

Var[Xt]  -  p2t  Var [Xq]  +  X6*-  (l-p2t)  • 

Thus  X  is  indeed  nonstationary.  However,  if  one  assumes  that  Xq  has  finite 
mean  and  variance  and  then  lets  t  go  to  infinity,  we  have 
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lim  E[Xt]  -  Xy  v 

t— >00  y 

lim  Cov[Xt,  X  ]  -  ^  e's/Y 
t— >00  *’  *’  " 


(2.1) 


lim  Var  [X  ]  -  X6£ 

t~>CD  ‘ 


(2.2) 


and  X  is  asymptotically  covariance  stationary  (see  Parzen  (1962a)).  Note 
that  letting  t  go  to  infinity  is  merely  saying  that  we  are  observing  the 
process  after  it  has  been  going  on  for  a  long  time. 

3.  CONDITIONAL  LEAST  SQUARES  ESTIMATORS  OF  DECAY  RATE  AND  AVERAGE  HEIGHT 
Since  Xt  is  an  asymptototically  stationary  autoregressive  process  of 

A  A  A  A 

order  one,  we  can  obtain  consistent  estimates  p,  y  *  -l/log(p),  and  o ~ 

A 

of  p,  y,  and  a2  =  Var[etl  by  choosing  p  as  the  value  of  p  minimizing 


s(P)  =  l  l(xt+1  -  X)  -  p(Xt  -  X)]: 
t=l 


(3.1) 


and  letting  o2  =  S(p)/(T-1) 

Unfortunately  this  procedure  gives  no  information  about  A  and 

However,  it  does  motivate  the  use  of  conditional  least  squares. 

If  {Zt,  t=0 ,+l , . . . }  is  an  ordinary  Gaussian  autoregressive  process  of 

order  one,  ijje.  Zt  =  a  +  et  where  E(et)  -  0  and  Cov(et,es)  -  6t_go2 

for  the  Kronecker  delta  function  6v»  then  E[Zt|  Zt_^,...,Z^]  *  a  • 

T-l 

Thus  the  minimization  of  ]»  {  Z^  -  a  Zt)2  can  be  thought  of  as  minimizing 

the  sum  of  squares  of  zt+^  from  its  conditional  expectation  given  all  previous 
Z's.  If  Z  is  not  Gaussian,  the  procedure  provides  the  coefficients  of  the 
best  linear  approximation  to  this  conditional  expectation.  In  our  situation 
we  need  to  actually  incorporate  the  conditional  expectation  into  the  function 
to  be  minimized  so  that  estimators  of  X  and  can  be  obtained.  Thus  we 
next  derive  this  conditional  expectation. 


it  is  Markov  and 


Since  X£+^  is  only  a  function  of  and  et+^, 
E[Xt+i!  Xx . Xt]  -  E[Xt+1|  Xt].  Further 

E(Xt+1|Xt1  -  E[pXt  +  «m|  Xt] 

-  PX,  +  Eltt+1  |Xt)  . 


But  Is  independent  of  Xt  so  that 

Elet+liXt3  *  E[et+13  "  xV(l-e'1/Y)  . 

Thus  the  conditional  least  squares  estimators  of  A,y^,  and  y  are  chosen 
to  minimize 


T-l 

1 

t-l 


J1  {Xfc+1  -  e‘1/Y  Xt  -AuY(l-e_1/Y)}2  . 


Unfortunately  A  and  yy  only  enter  this  function  in  their  product  and  thus 
they  are  not  identifiable.  However,  we  do  have  another  estimator  of  A  if 
we  can  count  how  many  episodes  occur  in  our  data.  Thus  we  will  assume  this 

a 

is  true  and  let  A  -  K/T  where  K  is  the  number  of  episodes  in  data  X^,...,^. 

A  A 

Then  we  let  B  -  Ay^  and  find  y  and  B  to  minimize 

T_i 

S(y,B)  -  l  {X  _  -  e-1/Y  X  -  By(l-e"1/Y)}2  .  (3.2) 

t-l 

To  find  the  properties  of  these  estimators  we  have  the  following  Theorem, 
adapted  in  a  straightforward  way  from  Theorem  3.2  of  Klimko  and  Nelson  (1978) 


Theorem 


The  estimators  y,  B  obtained  by  minimizing  (  3.2)  are  consistent  and 

[> 


•*(•;'  g)-5*  V?-0’  v'1> 


where  a2  -  Varte  1  -  Afcf  (1-e  2^Y) ,  and 
e  t  i 


1—  ,  1.  ■  .IL. 


a2  +02)  +  2ab 6y  +  b2  ac6y  +  be 

e2 

ac6y  +  be 


» 


with 

a  =  y"2e"1/Y  ,  b  =  B[l-Y_1  e_1/Y(l+y)]  , 

c  =  y(l-e  ly,Y)  . 


Proof  Theorem  3.2  in  Klimko  and  Nelson  (197-8)  shows  that  under  certain 


regularity  conditions  that  if  a  stochastic  process  Xfc  is  Markov  and 

asymptotically  stationary  and  g(0)  =  E£xt+i  |  Xt]  is  a  function  of  parameters 

T  A 

0  *  (0^,  .  ..,0^)  then  the  conditional  least  squares  estimators  0  of  0  ,  basei 


on  observations  X^,...,Xr,  satisfy 


Sr  (©  -  ©)  — >  n  (o,  v-1wv  x) 


r  - 


where  the  (j,k)th  elements  of  W  and  V  are  given  by 


W.,  =  E 


k+i-Mxt+1|V,  ^|xttl|xtf 

L  3ek 


and 


VE 


3E[Xt+1|Xt]  3E[Xt+1|Xt] 


30 


j 


30. 


j ,k=l, . . . ,r.  The  unconditional  expectations  in  these  expressions  are  the 
asymptotic  expectations. 

In  our  situation. 


3E[Xt+1|Xt] 


aXt_i  +  b 


(3.3) 


3y 

3E[Xt+l  |Xt] 


36 


y(l-e  A/^) 


(3.4) 


where  a  ■  y  2e  and  b  «  8[1  -  y  ^e 


(1+y) ]  .  Further, 


(Xt+1  -  E(Xt+1|Xtir  -  (pxt  +  €t+1  -  pxt  -  By (1-e 


‘  !et+l  '  E‘et+1>>‘ 


which  is  independent  of  3E[Xt+1  Ix^/Sy  and  3E[Xfc+^  | X^  J /30 .  Then 

W  -  o2V  and  thus  V_1WV_1  =  a2V_1.  Now  from  (2.1),  (2.2),  (3.3),  and  (3.4), 
e  e 


V22  =  E[cz]  =  cz  , 

V12  =  V21  =  E[acX  ,  +  be] 


=  acSy  +  be 


Vn  =  a2E[X2_1]  +  2  ab  ElX^]  +  b5 


a^y  (y  +  B2)  +  2abSy  +  b^ 


since  E[X2_2 ]  =  Var  [X^]  +  (ElX^])2  -  6+  A2p2  y  =  y(^-+02). 


To  minimize  the  nonlinear  function  (  3.2)  we  can  use  initial  values 
Y*  from  (3.1  )  and  8*  =  Ap*  where  A  *  K/T  as  above  and  p*  is  the 

average  of  the  heighths  observed  in  the  time  interval  containing  the 
start  of  an  episode.  Thus  p^  will  be  biased  downward  but  8  shduld  suffice 
as  a  starting  value  for  6. 

To  obtain  asymptotically  correct  confidence  intervals  for  y  and  8  we 
can  evaluate  o2  V-  with  o2  from  (3.1  )  and  A,  p^,  y,  and 

a  a  n  / 

&  «  a2/(A^-(l-e  /v)  replacing  a2,  A,  p  i  -,  y»  and  5. 

To  obtain  confidence  intervals  for  p^  we  can  use  the  facts  that 
-P 

p^  *  8/ A,  that  A-*-A  and  that  0  is  asymptotically  normal.  Thus  p^  is 
asymptotically  normal  with  asymptotic  variance  1/A2  times  the  asymptotic 
variance  of  0. 
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We  illustrate  these  ideas  in  the  next  section. 


4.  A  WORKED  EXAMPLE 

In  table  1  we  give  the  values  of  LH  in  a  cow  at  10  minute  intervals  for 
24  hours.  We  note  the  occurrence  of  K=8  episodes  in  the  T=144  points  with 
the  episodes  occurring  during  the  13th,  42nd,  58th,  69th,  89th,  100th,  115th, 
and  132  nd  time  intervals.  Thus  we  estimate  X  by  X  m  K/T  ■  8/144  *  .056 
occurrences  per  10  minutes  or  8  per  day  and  using  the  Poissoness  of  K,  we 
get  an  exact  90%  confidence  interval  for  XT  as  the  interval  [4.0,  14.5]. 

The  assumption  of  Poisson  occurrences  can  also  be  tested  by  (see  Parzen 


(1962),  p.  141)  comparing  the  statistic 


Z  * 


KT 

2 


618  -  576 
117.58 


.36 


with  the  critical  values  of  the  N(0,  1).  Thus  the  LH  data  is  clearly 
consistent  with  the  hypothesis  of  Poisson  occurrences. 

The  initial  estimator  of  y  given  by  the  ordinary  autoregressive  analysis 


is  y 


2.50  which  is  consistent  with  the  observation  that  it  takes  about 


8  or  10  intervals  for  the  episode  to  decay  back  to  the  baseline.  The  initial 


estimator  p . 

y 


20.5  and  thus  8  =  Xp..  =  1.14.  Then  using  the  IMSL  Subroutine 


ZXSSQ  to  minimize  (3.2  ),  we  obtain  y  =  2.50  and  8  *  1.90  and  thus 
*  ** 


p  *  33.93,  a  considerable  increase  over  p  . 

y  y 

The  estimator  o2  of  a2  =  X<S  ^  (1-e  ^Y)  obtained  from  the  autoregressive 
e  e  2 

analysis  is  a2  *  29.57,  from  which  we  obtain  6  *=a2/[Xy/2(l-e  ^Y)]  ■  386.59. 

A  AAA  A  A  A 

Then  substituting  a2,  6,  X,  y,  and  8  into  o2  and  V,  we  have  that  y  and  8  are 


approximately  normal  with  means  y  and  8  and  variances  .00153  and  .00123 


respectively,  from  which  we  have  asymptotic  95Z  confidence  intervals 

(2.42,  2.58)  and  (1.83,1.97)  for  y  and  8  respectively.  Finally  the  asymptotic 

A.  A  A 

variance  of  is  approximately  Var  (8)/X2  * 
confidence  interval  for  p^  is  (32.67,  35.19). 


.39  and  thus  an  asymptotic  95Z 
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Table  1.  LH  Levels  Every  10  Minutes  for  24  Hours 


No. 

Val. 

No. 

Val. 

No. 

Val. 

1 

8 

25 

7 

49 

11 

2 

7 

26 

7 

50 

10 

3 

7 

27 

7 

51 

10 

4 

8 

28 

7 

52 

10 

5 

8 

29 

6 

53 

8 

6 

7 

30 

6 

54 

9 

7 

5 

31 

6 

55 

7 

8 

7 

32 

5 

56 

8 

9 

6 

33 

4 

57 

10 

10 

6 

34 

4 

58 

28 

11 

6 

35 

6 

59 

25 

12 

7 

36 

4 

60 

20 

13 

22 

37 

6 

61 

18 

14 

36 

38 

7 

62 

16 

15 

24 

39 

5 

63 

13 

16 

19 

40 

6 

64 

12 

17 

17 

41 

6 

65 

11 

18 

13 

42 

22 

66 

9 

19 

12 

43 

36 

67 

11 

20 

13 

44 

27 

68 

9 

21 

9 

45 

22 

69 

26 

22 

9 

46 

18 

70 

29 

23 

9 

47 

15 

71 

24 

24 

8 

48 

14 

72 

18 

No. 

Val. 

No. 

Val. 

No. 

Val 

73 

15 

97 

9 

121 

10 

74 

13 

98 

8 

122 

11 

75 

12 

99 

11 

123 

10 

76 

10 

100 

31 

124 

8 

77 

9 

101 

26 

125 

8 

78 

9 

102 

21 

126 

8 

79 

9 

103 

18 

127 

9 

80 

9 

104 

16 

128 

8 

81 

8 

105 

13 

129 

8 

82 

7 

106 

12 

130 

8 

83 

7 

107 

11 

131 

8 

84 

8 

108 

9 

132 

14 

85 

7 

109 

10 

133 

25 

86 

7 

110 

9 

134 

20 

87 

8 

111 

9 

135 

16 

88 

6 

112 

7 

136 

15 

89 

42 

113 

9 

137 

11 

90 

26 

114 

9 

138 

11 

91 

21 

115 

32 

139 

10 

92 

18 

116 

25 

140 

10 

93 

19 

117 

21 

141 

13 

94 

13 

118 

22 

142 

11 

95 

11 

119 

18 

143 

11 

96 

12 

120 

13 

144 

10 

